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Abstract In this paper we present a simple and robust method 
for self-correction of camera distortion using single images 
of scenes which contain straight lines. Since the most com- 
mon distortion can be modelled as radial distortion, we il- 
lustrate the method using the Harris radial distortion model, 
but the method is applicable to any distortion model. The 
method is based on transforming the edgels of the distorted 
image to a 1-D angular Hough space, and optimizing the dis- 
tortion correction parameters which minimize the entropy 
of the corresponding normalized histogram. Properly cor- 
rected imagery will have fewer curved lines, and therefore 
less spread in Hough space. Since the method does not rely 
on any image structure beyond the existence of edgels shar- 
ing some common orientations and does not use edge fitting, 
it is applicable to a wide variety of image types. For instance, 
it can be applied equally well to images of texture with weak 
but dominant orientations, or images with strong vanishing 
points. Finally, the method is performed on both synthetic 
and real data revealing that it is particularly robust to noise. 



Keywords Radial distortion • Camera distortion • Plumb- 
line constraint 



1 Introduction 

1 . 1 Camera Calibration and Distortion Correction 

Camera calibration addresses the problem of finding the pa- 
rameters necessary to describe the mapping between 3-D 
world coordinates and 2-D image coordinates. This can be 
divided into the determination of extrinsic and intrinsic pa- 
rameters. The extrinsic parameters are necessary to relate 
an arbitrary 3-D world coordinate system to the internal 3- 
D camera coordinate system, where the z-axis is typically 
taken to be along the optical axis. The intrinsic parameters 
are independent of the camera position, and typically in- 
clude a focal length, an offset relating pixel coordinates cen- 
ter to the optical center, a scale factor describing the pixel 
aspect ratio, and one or more parameters describing non- 
linear radial and tangential distortions. Detailed examples 
of this process can be found in 11311161 , among others. 

In this work we address the determination of the non- 
linear distortion parameters, with the assumption that the 
other intrinsic parameters, usually addressed assuming a pin- 
hole camera model, can be found after the initial non-linear 
distortion correction. In order to do this a model must be 
selected, as well as a method of estimating the model pa- 
rameters. 
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1 .2 Summary of Related Work 

7.2.7 Distortion Correction Models 

A general equation for distortion correction is given by 

/i(xO=/o(x), x' = D(x), (1) 

where h is the output (corrected) image, /q is the input (dis- 
torted) image, and x is a pixel location in the input image 
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that is mapped to location x' in the output image by distor- 
tion correction function D(x) : ^ R^. 

If the distortion is assumed to be isotropic, and strictly 
radial, then D(x) can be written as 



D(x)=/(p)f + c, 



(2) 



where c is the center or radial distortion, the radius vector 
r = [ri r2]^ = X — c, the radius p = ||r||, the normalized ra- 
dius vector r = r/p and / is a scalar function of p. 

Several models for distortion correction have been used 
previously. In terms of /, the models include the polynomial 
model |[T3ll2l, where 



/(p) =p{\^hp^^k2P^^hp^ + ...), 

and the Harris model |27|, where 



(3) 



(4) 



More complex models such as fT4l include additional 



tangential distortion 



D(x) = 



7^1 (P^ + 
2pirir2 



2r\)^2p2rir2 



-P2{P 



-/(P)r + c, 



(5) 



to account for decentering of the optical system (where p\ 
and p2 are the parameters). More recently the rational func- 
tion model, where polynomial and perspective transforms 
are combined into a rational polynomial | 5 , 19 | has been pro- 
posed. 

All of these models can have adequate performance de- 
pending on the characteristics of the camera optics. In the 
following we limit our discussion and implementation to 
a single model in order to bound the scope of the paper, 
though it should be noted that the method we present is ex- 
tensible to any model. The Harris model is chosen for its 
reasonable performance, single parameter, and ease of in- 
version (achieved by simply negating 7). 

7.2.2 Parameter Estimation Methods 



from multiple positions or views of a particular scene. This 
is problematic in instances where access is no longer avail- 
able to the camera, but only to the resulting imagery. 

As a result, a number of methods have been developed 
which can operate on single views but which make use of 
common structure in images such as vanishing points [31, 
higher-order correlations in the frequency domain ifTTIl and 
straight lines |4|. 

The methods relying on the existence of straight lines, 
the "plumb-line" methods, were pioneered by Brown 141 . 
Brown used a test-field with actual strung plumb-lines, but 
in general these methods do not require knowledge of the 
locations of the lines. In particular, white plumb-lines were 
strung across a black background, with the plumb-bobs im- 
mersed in oil for stability. The photographic plate was ex- 
posed twice with the camera rotated about the optic axis by 
90°, giving a nominally square grid of lines. A number of 
points along the lines recorded on the photographic plate 
were measured using a microscope (a Mann comparator). 
Calibration was then performed by minimizing the least- 
squares error in the image between distorted straight lines 
and the measures points on the photographed plumb-lines. 

Plumb-line methods are applicable to many image types, 
because nearly all scenes containing man-made structures 
have a large number of straight lines. Plumb-line methods 
generally rely on the process of optimizing the distortion 
correction parameters to make lines that are curved by ra- 
dial distortion straight in the corrected imagery. The lines 
can be manually selected, as in |26|, or they can be found 
automatically, as in 11711251161 by detecting edgels and linking 
them in to line segments. 

The objective function for optimization can be formu- 
lating by undistorting the line segments and measuring the 
straightness by fitting a straight line 1261171. Alternatively, 
the distortion model can be chosen so that straight lines be- 
come specific family of curves, such as conies |6| or cir- 
cles I25I . The distortion can then be found by fitting these 
curves to the distorted line segments. 



A variety of methods exist for estimating camera distortion 
correction model parameters. Earlier efforts relied on im- 
agery with artificially created structure, either in the form of 
a test-field, populated with objects having known 3-D world 
coordinates, or using square calibration grids with lines at 
constant intervals I13[|16[ |2l . Alternative approaches do not 
require artificially created structure, but used multiple views 
of the same scene. The calibration technique makes use of 
constraints due to known camera motion (for instance rota- 
tion) [23 1, known scene geometry such as planar scenes |21 1 
or general motion and geometry constrained with the epipo- 
lar constraint ll24l[ni5il. 

These approaches required access to the camera in or- 
der to perform a specific operation, such as acquiring views 



2 Algorithm Overview 

We propose a method that is simple and robust to high levels 
of noise, as shown in the results section. In our algorithm we 
calculate all image edgels, and then transform these into a 
one-dimensional Hough space representation of angle. This 
creates an orientation histogram of the edgel angles. In this 
form, curved lines will be represented at a variety of angles, 
while straight lines will be found only at one. Therefore, we 
optimize the model distortion parameters which minimize 
the entropy (or spread) of the Hough space angular repre- 
sentation. The individual steps are: 

1. Find salient image edgels, with normal vectors. 
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(a) (b) (c) 

Fig. 1 Salient edge detection using tensor voting, (a) The original im- 
age, (b) the gradient magnitude, (c) the edge saliency. 



we split the image up in to a regular grid, and perform sub- 
sampling independently within each grid cell, so that each 
cell contains the same number of edgels. 

2.2 Distortion correction model 
2.2.7 Anisotropic extension 



2. Perform a distortion correcting transformation to the edgels. 

3. Compute the 1-D angular Hough transform. 

4. Compute an objective function defined as the spread (en- 
tropy) of the 1-D Hough transform. 

5. Optimize the transformation parameters to minimize the 
entropy/spread based objective function, iterating steps 
2-4, 

6. Use the optimized transform parameters to map the input 
image to a corrected output image. 

Note that step 1, finding the edgels, is only required 
once, due to the fact that the edgels, rather than the underly- 
ing image, can be transformed. This, and the other steps of 
the process, are described in further detail in the following 
respectively enumerated subsections. 

2.1 Salient edge extraction 

The structure that we are making use of in this paper con- 
sists of long, straight edges, which appear as long, smoothly 
curved edges in the input image. 

Simply performing an edge detection does not result in 
these features being dominant, so to enhance the long edges, 
we use tensor voting 1 15 1 on the dense gradient image. The 
gradient is produced by finite differences, and voting is made 
in proportion to the gradient magnitude. We use the stan- 
dard kernel for smooth, circular curves, specifically the ker- 
nel given by the equations on page 60 of ifTSl . 

For edges, each pixel is represented by a 2 x 2 positive 
semi-definite matrix with eigenvalues Ai and A2. We wish to 
find points which are sufficiently 'edgey' so we use the edge 
saliency function 0, given by: 

= maxAi, A2 — ^minAi, A2. (6) 

This detects points where the edge in the main orientation is 
significantly brighter than other edges. We have found that 
e = 2 produces good results. An example of these are shown 
in Figure[T] The normal of these edgels is given by the eigen- 
vector corresponding to the largest eigenvalue. 

The final stage is to discard the non- salient edges by 
thresholding at = 0. For computational efficiency of the 
later stages, the set of edges passing the threshold is subsam- 
pled. Typically, we keep 100,000 edgels. In order to prevent 
the result being dominated by one region with strong edges. 



As mentioned previously, for the basic radial distortion cor- 
rection model, we have based our distortion correction model 
on the Harris model defined in Equation]?] In order to extend 
the model's flexibility, we add an anisotropic component, so 
that the model becomes: 

D(x) = f/(p)(l+^(0))+c, (7) 

where /(.) and g{.) are the distortion functions. 

In particular, we define the anisotropy function, g{6) to 

be: 

g{e) = aisin(e + V^i)+^/2sin^(0 + V^2) + .-- 

= Z7isine + Z72COS0 + (/73sine + Z74COS0)^ + ..., (8) 

which can be rearranged as a Fourier series. This general- 
purpose formulation has no discontinuity as the angle wraps 
around and conveniently avoids the use of trigonometric func- 
tions, since cos 6 = ri/p, etc. 

2.2.2 Edgel transformation 

In order to evaluate the cost function it is necessary to de- 
termine the edgels' orientation, as discussed in the previous 
section. If it were necessary to regenerate these by comput- 
ing a completely new image for each new set of distortion 
model parameters, however, the optimization process (in- 
volving numerous evaluations of different sets of model pa- 
rameters) would be prohibitively slow. Furthermore, errors 
due to resampling would be introduced when the image was 
undistorted at each iteration. 

Fortunately, the edgels need to be determined only once, 
from the input image, due to the fact that the distortion cor- 
recting transformations can be applied directly to the edgels, 
provided that the Jacobians of the transformations are known. 
In order to see this, consider an edgel in the input image as 
part of a parameterised curve l(^) in R^. The distortion cor- 
rection transformation, D, can be applied to this space to 
create a new line: 

m(0=D(l(0). (9) 
The line tangents are then found by differentiating to be: 
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Fig. 2 Images of curved and straight line images and tlie correspond- 
ing 1-D Hough transforms. 



where the Jacobian, J, is given by: 



dli dh 
dD2 dD2 
dli dh 



(11) 
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The edgels include normals n defined at discrete points. The 
transformed edgel normals, h, can then be found by rotating 
the normals by 90° (making them tangents), transforming 
the tangents by J and rotating the transformed tangents back 
to be normals: 



(12) 



Note that we do not parameterise the line with a func- 
tion. The line and its normal is known (and used) only at a 
discrete set of points, specifically where the edgels are de- 
tected. This means that l(^) and n(^) can be evaluated at ev- 
ery value of t we require. Since the edgel detection process 
also provides the normals, J is only a function of the distor- 
tion model, and is therefore computed analytically from the 
definition of D. The derivation of J for the Harris model is 
given in Appendix [A| 



2.3 1-D Angular Hough transform 

The Hough transform is a technique for finding lines in im- 
ages |9l. It is sufficiently well-known so that only a trun- 
cated explanation will be given here. The basis of the tech- 
nique is the transform of a line to a point in "Hough" space. 

By way of example, if a line is defined by y = mx^h, 
then it can be represented by a single point in a 2-D Hough 
space of [m] x [h] . An edgel, which contains both a normal 
vector (and thus the slope of a line) and a discrete point 
(which a corresponding line would pass through), can also 



be mapped to a single point in Hough space. In practice the 
mx + b formulation is unwieldy, so 



p = XCOS0 + jsinO 



(13) 



is used, where p is the shortest distance from the line to the 
origin, and 6 is defined by the vector normal to the line. 

The Hough space can be quantized into discrete bins, 
and each edgel can be assigned to a particular bin. Given a 
set of edgels, support for the existence of a particular line 
is then indicated by the accumulation of a large number of 
edgels in the corresponding bin. The edgels from a curve 
end up in adjacent bins, resulting in a diffuse cluster of non- 
zero bins, while edgels that are all on the same line end up 
in exactly the same bin. This important fact motivates the 
definition of our objective function, as explained in the next 
section. 

In general we anticipate the presence of a number of 
parallel lines in the scene, and would like to utilize their 
reinforcement of each other. We therefore marginalize the 
[p] X [0] Hough space into a 1-D space by summing over 
the p values for each B. An example of some images with 
curved lines and straight lines, and their corresponding 1-D 
Hough transforms, are shown in Fig. [2] Note that the two 
principal line directions are mapped to two corresponding 
peaks in Hough space. 



2.4 Entropy-based objective function 

The radial distortion correction method presented here is 
motivated by the observation that curved lines map to spread 
out peaks in Hough space, while straight lines map to a sin- 
gle bin. Therefore, it is desirable to have an objective func- 
tion that measures this spread. In information theory this 
quality is represented by entropy |22|. We have therefore 
normalized the 1-D Hough representation, and treat it as a 
probability distribution. The objective function is then: 



C(H) = -£p(//,)/og2(p(H,)), 

b=l 



(14) 



where H is the discretized 1-D angular Hough transform of 
the edgels of an input image for a given set of model pa- 
rameters, p{Hiy) is the value of a given normalized Hough 
bin, and B is the number of bins. Minimizing C therefore 
minimizes the spread. 



2.5 Optimization 

The cost function has many local minima, so to optimize 
it effectively a reasonable strategy is needed. As a broad 
overview, we have found that the MCDH (Monte-Carlo down- 
hill) strategy works best. Essentially, starting parameters are 



5 



selected at random, and the downhill simplex algorithm flTll 
(specificially the variation given in 1 18 1) is used to optimize 
the cost from the selected starting points. The best result 
is then selected. The downhill simplex algorithm is iterated 
until the residual drops below a threshold (10~^^), or 1000 
iterations occur, whichever is sooner. The residual is com- 
puted as: 

1 1 Sbest ~ Sworst 1 1 2 



Uncorrelated noise 



||Sbest||2 

where Sbest is the location of the best (smallest) point on the 
simplex, and Sworst is the location of the worst (largest point) 
on the simplex. 

In order to implement this method effectively, several 
techniques are needed. Parameters computed in the Hough 
transform will not lie exactly at the bin centres. If the param- 
eter is simply placed in the closest bin, then the quantization 
caused by this causes flat areas to exist in the cost function. 
The flat areas can easily confound the optimizer, so instead, 
the neighbouring bins are incremented using linear interpo- 
lation. 

The optimizer is most effective when the parameters be- 
ing optimized have similar orders of magnitude, since other- 
wise roundoff error will cause errors in the computation of 
new simplexes. In this case, a simple scale based on a priori 
knowledge of approximate parameter values suffices. 

For radial distortion, c is roughly in the middle of an 
image and is of the order 10^ to 10^ pixels. Given Equa- 
tion [4j it is reasonable to expect 7 to be approximately in the 
range ±Pmax' where Pmax ^ |c| and c is the centre of the im- 
age. Therefore, in the optimization, we solve for j3, where 
j3 = lOOyp^jjx' which brings j3 in to the same range as c. 

Similarly, we have found that the values of the anisotropic 
coefficients, are of order 10~^, so we instead solve 

for 4, where di = 10^ bf. The resulting optimization is per- 
formed over the vector: 



[ci,C2,P,di,d2,d3,d4,d5,i 



(15) 



The random selection of the starting parameters must be 
based around some knowledge of the values. The centre of 
radial distortion is usually near to the centre of the image, so 
we draw the initial centres from a normal distribution cen- 

Icl 

tred at c with a standard deviation of ^ . Given the approxi- 
mate values of 7, and therefore j3 , we draw initial values of 
P from ^(0, 100). The initial anisotropy parameters are set 
to zero. For our application and these parameters, we have 
found that 120 MCDH iterations is sufficient. 



2.6 Generating the output image 

In order to generate the output (corrected) image the follow- 
ing procedure is needed: 
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Fig. 3 (Left) Illustration of the simulated data (shown with 100 points 
per line). (Right) Performance of the algorithm on the simulated data, 
with Y= 10~^ and 7 = 2 x 10~^. The median value of gamma is shown 
with error bars at the 10* and 90* percentile. 



1. For each point, x', in a grid defined on the output image 
/i , find the corresponding point in the distorted input im- 
age Iq given by (x'). 

2. Copy the pixel back, performing the assignment: 

/i(x'):=/o(D-i(x')). 

Since images are discrete, interpolation will be required 
to find /i(D~^(x')). Bicubic interpolation has been used to 
generate the results shown in the following section. 



3 Results 

3.1 Synthetic tests 

In order to compute the sensitivity to various kinds of noise, 
we create a number of synthetic images, using known radial 
distortion levels, and compute the resulting radial distortion 
correction parameters using our technique. We investigate 
the two correction models, namely the strictly radial cor- 
rection of the unmodified Harris model, and the anisotropic 
extension discussed previously. The images are designed to 
test the performance of the algorithm in the presence of mea- 
surement error in the plumb-lines, and both uncorrelated and 
correlated non-plumb-line data. 

The 250 X 250 pixel test images consist of five good 
quality lines (no point can be within 60 pixels of the centre 
of the image), consisting of 10 points on each line, drawn 
with barrel distortion and with random orientation errors. 
Then (in the case of uncorrelated noise), a number of ran- 
domly placed, randomly oriented points are added. In the 
case of correlated noise, a number of randomly placed, sized 
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and oriented ellipses (with the same approximate point spac- 
ing as the lines) are added. Our algorithm is then applied to 
the test image to estimate the parameters. 

For every selected noise proportion for the simulation, 
200 test images are created for each of two different values 
of 7. The results of these, and an illustration of the simulated 
data is shown in Figure |3] As can be seen, the technique is 
very robust to noise, producing good results with up to 70% 
contamination and interestingly it does not perform much 
worse in the case of correlated noise. Furthermore, although 
the anisotropic extension (in this case) introduces six extra 
parameters to the optimization, this has no significant effect 
on the stability of the technique in high noise situations. 

3.2 Example images 

The results on two real images are shown in Figure]?] In the 
corrected images, real- world straight lines have been anno- 
tated with lines to show that there is no significant amount 
of curvature. The upper image shows the technique operat- 
ing on an image where the lines belong to vanishing points 
(note that some distortion was artificially added to this im- 
age in addition to the camera distortion present). Note the 
presence of strong curved edges along the boundaries of the 
cars and parts of the building. 

The lower image shows the results on an aerial image 
of a city. Although there are two principal directions in this 
image, the contrast on the principal edges is low compared 
to many of the other features. Additionally, the anisotropy 
extension was required to correct the very strong distortion 
in the lower left comer. 



3.3 Comparison to other techniques 

To gauge the accuracy of our method, we have compared it 
to a technique based on structured scenes. In particular we 
have used the technique of 1 10 |, which uses multiple images 
of a planar grid of squares to determine the intrinsic and 
extrinsic camera calibration parameters. This system is able 
to calibrate cubic, quintic, Harris and Harris with unit aspect 
ratio camera models. We use the last listed model, since this 
matches the camera model in the paper. In particular, we use 
this program to optimize the model: 
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where (u^v) is the optic axis, is a coordinate in the 

image, and (x^, jc) is a coordinate in ideal, normalized cam- 
era's image plane. Since the technique is a structured scene 
technique, it also estimates the focal length of the camera, 
/. To translate from one model to the other, set 7 = ^ . Note 



that this assumes that the centre of radial distortion is at the 
optic axis of the camera, which is a reasonable assumption. 

In order to measure the quality of the camera calibra- 
tion, we localise a known 3D model in an image, and mea- 
sure the root mean square (RMS) error between the control 
points on the model and the measured points in the image. 
The model consists of 13 by 9 black and white squares in a 
chequerboard pattern on a flat plane. The model is warped 
using a homography (the parameters of which must be de- 
termined), then by the calibrated radial distortion function 
and rendered in to the image. We then search normal to the 
rendered edges, looking for edges in the image. We define 
the error to be the distance offset in pixels between the ren- 
dered model edge and the measured edge position in the im- 
age The parameters of the homography are then adjusted 
using the downhill- simplex algorithm to minimize the RMS 
error. In particular, we search 15 pixels in either direction, 
and take the point with the highest gradient, measured us- 
ing the kernel [—1 — 2 — 10121]. If the highest gradient 
does not exceed some threshold, or is not a local maxima, 
then the edge search is discarded and does not contribute to 
the RMS error. Quadratic interpolation is then used to find 
the edge position to subpixel accuracy. The search distance 
is sufficiently small that data association is trivial, and ro- 
bust fitting is not required. A more detailed treatment of very 
closely related systems for aligning models is given in (121 
8|. 

The results of this test are shown in Figure [5] As ex- 
pected, the technique which uses many (in this case approxi- 
mately 100) images of a structured scene performs best. The 
performance of our method varies somewhat with the con- 
tents of the scene, but performs well on certain images. One 
interesting property of our technique is well illustrated. The 
property is that edgels from multiple, unrelated images can 
be aggregated to provide more straight edges. In the case 
shown, the resulting calibration is superior to the calibra- 
tion generated from a single image of a carefully constructed 
grid. 



4 Discussion and conclusions 

In this paper, we have presented a new, simple and robust 
method for determining the radial distortion of an image us- 
ing the plumb-line constraint. The technique works by first 
extracting salient edgels and then minimizing the spread of 
a ID angular Hough transform of these edgels. The tech- 
nique is simple and because no edge fitting is performed, the 
technique is very robust to the presence of noise. Further- 
more, the technique is more generally applicable than other 
plumb-line techniques in that the lines used do not need to 
be continuous. The technique works on textures with prin- 
cipal directions, as illustrated by the aerial image of a city. 



Fig. 4 (Left) Original, distorted images. (Right) Images undistorted using our technique. Some lines which are straight in the world have been 
annotated with straight lines in the undistorted image. The images are of a building, showing strong vanishing points and an aerial image of a city. 



where the salient edge detection results in a large number of 
relatively small edge fragments. 

The proposed algorithm has a number of parameters: the 
parameters of the tensor voting kernel, the number of bins 
and the parameters of the optimization. In practice, the se- 
lection of these parameters are not critical, and indeed the 
same set of parameters was used for the simulated data, the 
example images and the test images shown. 

Our method is flexible in that it does not impose con- 
straints beyond the presence of one or more straight edges: 
it is not a requirement that the edges share vanishing points, 
or structure of any particular kind. It is not even a require- 
ment that the edgels belong to a related set of images. The 
technique can be equally applied to edgels from multiple 
images of unrelated scenes taken with the same camera pa- 
rameters. Finally, our method is widely applicable because 
it is, in terms of RMS error, able to produce a calibration 
to within three percentage points of a technique requiring 
access to the camera and structured scenes. 
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Fig. 5 Calibration results for an alternative method and some test images, including a calibration grid and some images of urban scenes. All images 
are 2048 x 1536 pixels. % In the uncalibrated case, the errors are sufficiently that the search radius had to be increased to 100 pixels. This results 
in the RMS error being an underestimate, since for large errors, the system sometimes finds a closer edge (so the measured error is too small), or 
fails to find the edge at all (so large errors are not included), t In some cases, the errors caused a large number of the correspondences to be missed 
with a search radius of 15 pixels, so it was increased to 36 pixels. For the last row of the table, the edgels from the two images were aggregated 
and treated as a single image. In our model, c = [c\ cj\ and for the libCVD model, u = c\ and v = C2. 
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A Derivation of J for the Harris model 

^1 

Taking f to be a column vector (i.e. f = ^ ): 

\r2_ 



-f/(p) 



^g(9(x)) 



Where: 



5f 

dm 

dp 
dp 
dx 



-dn_ 




1 




-hh 


dxi 
_ dx\ 




dh 
dx2 . 


~ P 







and: 



dx 



[ dxi dx2 J 



\dg_ dp_'] 
[ dxi dx2 J 



n r2 , 



The function is a function of sin 6 and cos 6, where: 

cos Q = r\/ p = r\, and 
sin = r2 /p = r2 , 

we can rewrite g{0(x)) as: 

g(f) =bxr2^ b2r\ + {hh + b^h f + {b^h + behf , 
so, 

dg, ^ dgjr) dr _ idg^^^ ^^(j.)] dr 
dx dr dx 

where: 

dg{r 



[ dh dh J dx 



dr\ 
dgir) 
dr2 



bx+2bi{bifx +bir2)+ibs{bsh +b(,h) 



-- b2 + 2b4(bif\ +bj,f2) + 'ib(,{bifi+b(,f2) 



(17) 

(18) 

(19) 
(20) 

(21) 



(22) 
(23) 



(24) 
(25) 

(26) 
(27) 



Substituting Equations [T8l[T9l|20l|25]|26l and[27lin to EquationfTTlgives 
the analytic solution for the J. 



